##################################################################################################
# "Natural Disasters, 'Partisan Retrospection,' and U.S. Presidential Elections" #################
# Boris Heersink, Jeffrey Jenkins, Michael Olson, and Brenton Peterson ###########################
# Panel Analysis #################################################################################
##################################################################################################

# load data ######################################################################################

# gasper and reeves data

  gr <- read.csv("gr_uncorrected.csv")
  
  gr_corrected <- read.csv("gr_corrected.csv")

# summary statistics ##############################################################################

  gr_correct_forsumstats <- gr_corrected[,c("inc.twpct","propertypercapl.6mo","disdecs.all.6mo","turndowns.6mo",
                                            "copart","swing","medincK","pres.vote.lag1","pres.vote.lag2")]
  
  stargazer(gr_correct_forsumstats,
            summary=T,
            covariate.labels = c("Incumbent Party Two-Party Vote Share",
                                 "ln(Damage per 10,000)",
                                 "Disaster Declarations, 6 Months",
                                 "Turndowns, 6 Months",
                                 "Co-Partisan County",
                                 "Swing County",
                                 "Median Household Income",
                                 "Incumbent Party Two-Party Vote Share, t-1",
                                 "Incumbent Party Two-Party Vote Share, t-2"),
            summary.stat=c("n","mean","median","sd","min","max"),
            out="./results/gr_sumstats.tex",
            title="Summary Statistics, Presidential Elections and Disasters, 1972-2004",
            label="gr_sumstats",
            font.size = "footnotesize")

# estimate models #################################################################################

# new model specification - replicate gasper and reeves

  gr_base <- felm(inc.twpct ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo +  medincK| countyfips + year | 0 | countyfips ,data = gr_corrected)
  
  gr_base_int_new <- felm(inc.twpct ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo*(copart+swing) +  medincK| countyfips + year | 0 | countyfips ,data = gr_corrected)
  
  gr_aug <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo +  medincK| countyfips + year| 0 | countyfips ,data = gr_corrected)
  
  gr_aug_int_new <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo)*(copart+swing) +  medincK| countyfips + year| 0 | countyfips ,data = gr_corrected)
  
  stargazer(gr_base,gr_base_int_new,gr_aug,gr_aug_int_new,
            dep.var.labels = c("Incumbent Party Two-Party Vote Share"),
            title="Table 2: Impact of Disaster Damage on Incumbent Vote Share, 1972-2004",
            label="gr_new",
            order=c(3,6,7,9,11,13,10,12,14,4,5,8,1,2),
            covariate.labels = c("Damage",
                                 "Declarations, 6 Months",
                                 "Turndowns, 6 Months",
                                 "Co-Partisan X Damage",
                                 "Co-Partisan X Declarations",
                                 "Co-Partisan X Turndowns",
                                 "Swing X Damage",
                                 "Swing X Declarations",
                                 "Swing X Turndowns",
                                 "Co-Partisan",
                                 "Swing",
                                 "Median Income",
                                 "Incumbent Party Vote Share, Lag 1",
                                 "Incumbent Party Vote Share, Lag 2"
            ),
            type="html",
            keep.stat = c("n","rsq"),no.space = T,
            star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
            notes.append = FALSE,notes.label = "",column.sep.width="0pt",
            table.layout ="-ld-#-t-as-n",
            out="./results/gasperreeves_new.html",
            notes="Note: Entries are linear regression coefficients with county-clustered standard errors in parentheses. All models include county and year fixed effects. <sup>*</sup>p<0.10, <sup>**</sup><0.05 (two-tailed test).")
  
  
# make figure with scaled variables
  
  gr_corrected_scaled <- gr_corrected
  gr_corrected_scaled[,c("inc.twpct","pres.vote.lag1","pres.vote.lag2","propertypercapl.6mo","medincK","disdecs.all.6mo","turndowns.6mo")] <- 
    scale(gr_corrected_scaled[,c("inc.twpct","pres.vote.lag1","pres.vote.lag2","propertypercapl.6mo","medincK","disdecs.all.6mo","turndowns.6mo")])
  
  gr_aug_int_new1 <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo)*(copart+swing) +  medincK| countyfips + year| 0 | countyfips ,data = gr_corrected_scaled)
  gr_aug_int_new2 <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo)*(contra+swing) +  medincK| countyfips + year| 0 | countyfips ,data = gr_corrected_scaled)
  gr_aug_int_new3 <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo)*(copart+contra) +  medincK| countyfips + year| 0 | countyfips ,data = gr_corrected_scaled)
  
  contra_fig <- as.data.frame(summary(gr_aug_int_new1)$coefficients[c("propertypercapl.6mo","disdecs.all.6mo","turndowns.6mo"),])
  contra_fig$var <- rownames(contra_fig)
  contra_fig$type <- "Contra-Partisan"
  
  copart_fig <- as.data.frame(summary(gr_aug_int_new2)$coefficients[c("propertypercapl.6mo","disdecs.all.6mo","turndowns.6mo"),])
  copart_fig$var <- rownames(copart_fig)
  copart_fig$type <- "Co-Partisan"
  
  swing_fig <- as.data.frame(summary(gr_aug_int_new3)$coefficients[c("propertypercapl.6mo","disdecs.all.6mo","turndowns.6mo"),])
  swing_fig$var <- rownames(swing_fig)
  swing_fig$type <- "Swing" 
  
  coef_fig <- as.data.frame(rbind(copart_fig,swing_fig,contra_fig))
  coef_fig$type <- factor(coef_fig$type,levels=c("Co-Partisan","Swing","Contra-Partisan"))
  coef_fig$var <- plyr::revalue(coef_fig$var,c("propertypercapl.6mo"="Damage",
                                               "disdecs.all.6mo"="Declarations",
                                               "turndowns.6mo"="Turndowns"))
  coef_fig$var <- factor(coef_fig$var,levels=c("Damage","Declarations","Turndowns"))
  
  cairo_ps("./results/fig3.eps",width=8,height=5,fallback_resolution = 1200)
  print(
    ggplot(data=coef_fig,aes(x=var,y=Estimate,ymin=Estimate-2*`Cluster s.e.`,ymax=Estimate+2*`Cluster s.e.`,group=type,colour=type))+
      geom_point(size=4,position=position_dodge(width=0.25))+
      geom_errorbar(size=2,width=0.25,position=position_dodge(width=0.25))+
      geom_hline(yintercept = 0,colour="indianred",size=2,linetype=2)+
      scale_colour_grey()+
      theme_minimal()+
      labs(colour = "Type of County")+
      xlab("Variable")+ylab("Standardized Coefficient")
  )
  dev.off()
  
  png("./results/fig3.png",width=8,height=5,units="in",res = 1200)
  print(
    ggplot(data=coef_fig,aes(x=var,y=Estimate,ymin=Estimate-2*`Cluster s.e.`,ymax=Estimate+2*`Cluster s.e.`,group=type,colour=type))+
      geom_point(size=4,position=position_dodge(width=0.25))+
      geom_errorbar(size=2,width=0.25,position=position_dodge(width=0.25))+
      geom_hline(yintercept = 0,colour="indianred",size=2,linetype=2)+
      scale_colour_grey()+
      theme_minimal()+
      labs(colour = "Type of County")+
      xlab("Variable")+ylab("Standardized Coefficient")
  )
  dev.off()

# original data models
  
  gr_base_orig <- felm(inc.twpct ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo +  medincK| countyfips + year | 0 | countyfips ,data = gr)
  
  gr_base_int_new_orig <- felm(inc.twpct ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo*(copart+swing) +  medincK| countyfips + year | 0 | countyfips ,data = gr)
  
  gr_aug_orig <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo +  medincK| countyfips + year| 0 | countyfips ,data = gr)
  
  gr_aug_int_new_orig <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo)*(copart+swing) +  medincK| countyfips + year| 0 | countyfips ,data = gr)
  
  orig_sg <- stargazer(gr_base_orig,gr_base_int_new_orig,gr_aug_orig,gr_aug_int_new_orig,
                       dep.var.labels = c("Incumbent Party Two-Party Vote Share"),
                       title="Impact of Disaster Damage on Incumbent Vote Share, 1972-2004: Original Data with Duplicates",
                       label="gr_orig",
                       order=c(3,6,7,9,11,13,10,12,14,4,5,8,1,2),
                       covariate.labels = c("Damage",
                                            "Declarations, 6 Months",
                                            "Turndowns, 6 Months",
                                            "Co-Partisan $\\times$ Damage",
                                            "Co-Partisan $\\times$ Declarations",
                                            "Co-Partisan $\\times$ Turndowns",
                                            "Swing $\\times$ Damage",
                                            "Swing $\\times$ Declarations",
                                            "Swing $\\times$ Turndowns",
                                            "Co-Partisan",
                                            "Swing",
                                            "Median Income",
                                            "Incumbent Party Vote Share, Lag 1",
                                            "Incumbent Party Vote Share, Lag 2"
                       ),
                       keep.stat = c("n","rsq"),no.space = T,
                       star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                       notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                       table.layout ="-ld-#-t-as-n",font.size="footnotesize",
                       notes="\\parbox[t]{0.75\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       county-clustered standard errors in parentheses. All models include county and year fixed effects. 
                       $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(orig_sg, sep = '\n', file = "./results/gr_corr.tex")

# five-bin measure of previous support  

  gr_corrected$quantile <- NA
  gr_corrected[gr_corrected$past3<35 & !is.na(gr_corrected$past3),"quantile"] <- 1
  gr_corrected[gr_corrected$past3>=35 & gr_corrected$past3<45 & !is.na(gr_corrected$past3),"quantile"] <- 2
  gr_corrected[gr_corrected$past3>=45 & gr_corrected$past3<=55 & !is.na(gr_corrected$past3),"quantile"] <- 3
  gr_corrected[gr_corrected$past3>55 & gr_corrected$past3<=65 & !is.na(gr_corrected$past3),"quantile"] <- 4
  gr_corrected[gr_corrected$past3>65 & !is.na(gr_corrected$past3),"quantile"] <- 5
  
  gr_base_int_quantile <- felm(inc.twpct ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo*(factor(quantile)) +  medincK| countyfips + year | 0 | countyfips ,data = gr_corrected)
  
  gr_aug_int_quantile <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (disdecs.all.6mo + turndowns.6mo+propertypercapl.6mo)*(factor(quantile)) +  medincK| countyfips + year| 0 | countyfips ,data = gr_corrected)
  
  gr_quantile_sg <- stargazer(gr_base_int_quantile,gr_aug_int_quantile,
                              dep.var.labels = c("Incumbent Party Two-Party Vote Share"),
                              title="Disaster Damage and Incumbent Vote Share, 1972-2004: Five-Bin Measure of Partisan Lean",
                              label="gr_quantile",
                              order=c(5,19,20,21,22,
                                      3,11,12,13,14,
                                      4,15,16,17,18,
                                      6,7,8,9,10,1,2),
                              covariate.labels=c("Damage",
                                                 "Damage $\\times$ 35$-$45",
                                                 "Damage $\\times$ 45$-$55",
                                                 "Damage $\\times$ 55$-$65",
                                                 "Damage $\\times$ $>$65",
                                                 "Declarations",
                                                 "Declarations $\\times$ 35$-$45",
                                                 "Declarations $\\times$ 45$-$55",
                                                 "Declarations $\\times$ 55$-$65",
                                                 "Declarations $\\times$ $>$65",
                                                 "Turndowns",
                                                 "Turndowns $\\times$ 35$-$45",
                                                 "Turndowns $\\times$ 45$-$55",
                                                 "Turndowns $\\times$ 55$-$65",
                                                 "Turndowns $\\times$ $>$65",
                                                 "35$-$45",
                                                 "45$-$55",
                                                 "55$-$65",
                                                 "$>$65",
                                                 "Median Income",
                                                 "Incumbent Party Vote Share, Lag 1",
                                                 "Incumbent Party Vote Share, Lag 2"
                              ),
                              keep.stat = c("n","rsq"),no.space = T,single.row = T,
                              star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                              notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                              table.layout ="-ld-#-t-as-n",#font.size="footnotesize",
                              notes="\\parbox[t]{0.85\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       county-clustered standard errors in parentheses. `Partisan Lean' is the mean incubment-party vote share from the prior
                       three presidential elections. All models include county and year fixed effects. 
                       $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(gr_quantile_sg, sep = '\n', file = "./results/gr_quantile.tex")
  
# continuous measure of prior support

  gr_base_int_cont <- felm(inc.twpct ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo*(past3) +  medincK| countyfips + year | 0 | countyfips ,data = gr_corrected)
  
  gr_aug_int_cont <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (disdecs.all.6mo + turndowns.6mo+propertypercapl.6mo)*(past3) +  medincK| countyfips + year| 0 | countyfips ,data = gr_corrected)
  
  gr_cont_sg <- stargazer(gr_base_int_cont,gr_aug_int_cont,
                          dep.var.labels = c("Inc. Party 2-Party Vote Share"),
                          title="Disaster Damage and Incumbent Vote Share, 1972-2004: Continuous Measure of Partisan Lean",
                          label="gr_cont",
                          order=c(5,3,4,10,8,9,6,7,1,2),
                          covariate.labels = c("Damage",
                                               "Declarations, 6 Months",
                                               "Turndowns, 6 Months",
                                               "Partisan Lean $\\times$ Damage",
                                               "Partisan Lean $\\times$ Declarations",
                                               "Partisan Lean $\\times$ Turndowns",
                                               "Partisan Lean",
                                               "Median Income",
                                               "Incumbent Party Vote Share, Lag 1",
                                               "Incumbent Party Vote Share, Lag 2"
                          ),
                          keep.stat = c("n","rsq"),no.space = T,
                          star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                          notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                          table.layout ="-ld-#-t-as-n",
                          notes="\\parbox[t]{0.7\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       county-clustered standard errors in parentheses. `Partisan Lean' is the mean incubment-party vote share from the prior
                       three presidential elections. All models include county and year fixed effects. 
                       $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(gr_cont_sg, sep = '\n', file = "./results/gr_cont.tex")
  

  cont_binning <- inter.binning(data=gr_corrected,Y="inc.twpct",D="propertypercapl.6mo",
                                X="past3",Z=c("medincK","disdecs.all.6mo", "turndowns.6mo",
                                              "pres.vote.lag1" , "pres.vote.lag2" ),
                                FE=c("countyfips","year"),full.moderate = T,na.rm=T,CI=F,cl="countyfips",
                                nbins = 3)

  pdf("./results/gr_cont_bins.pdf",width=8,height=5)
  print(
    
    plot.interflex(cont_binning,ylab="Marginal Effect of Disaster Damage",xlab="Incumbent Party Support")
    
  )
  dev.off()

# fully moderated models

  gr_corrected$newout1 <- gr_corrected$inc.twpct
  gr_corrected$newout2 <- gr_corrected$inc.twpct
  
  gr_corrected_scaled <- gr_corrected
  gr_corrected_scaled[,c("inc.twpct","newout1","newout2","pres.vote.lag1","pres.vote.lag2","propertypercapl.6mo","medincK","disdecs.all.6mo","turndowns.6mo")] <- 
    scale(gr_corrected_scaled[,c("inc.twpct","newout1","newout2","pres.vote.lag1","pres.vote.lag2","propertypercapl.6mo","medincK","disdecs.all.6mo","turndowns.6mo")])
  
  gr_aug_corr_copart <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + propertypercapl.6mo +  medincK +  disdecs.all.6mo + turndowns.6mo| countyfips + year| 0 | countyfips ,data = gr_corrected[gr_corrected$copart==1,])
  
  gr_aug_corr_swing <- felm(newout1 ~ pres.vote.lag1 + pres.vote.lag2 + propertypercapl.6mo +  medincK +  disdecs.all.6mo + turndowns.6mo| countyfips + year| 0 | countyfips ,data = gr_corrected[gr_corrected$swing==1,])
  
  gr_aug_corr_contra <- felm(newout2 ~ pres.vote.lag1 + pres.vote.lag2 + propertypercapl.6mo +  medincK +  disdecs.all.6mo + turndowns.6mo| countyfips + year| 0 | countyfips ,data = gr_corrected[gr_corrected$contra==1,])
  
  gr_aug_corr_copart_stand <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + propertypercapl.6mo +  medincK +  disdecs.all.6mo + turndowns.6mo| countyfips + year| 0 | 
                                     countyfips ,data = gr_corrected_scaled[gr_corrected_scaled$copart==1,])
  
  gr_aug_corr_swing_stand <- felm(newout1 ~ pres.vote.lag1 + pres.vote.lag2 + propertypercapl.6mo +  medincK +  disdecs.all.6mo + turndowns.6mo| countyfips + year| 0 | 
                                    countyfips ,data = gr_corrected_scaled[gr_corrected_scaled$swing==1,])
  
  gr_aug_corr_contra_stand <- felm(newout2 ~ pres.vote.lag1 + pres.vote.lag2 + propertypercapl.6mo +  medincK +  disdecs.all.6mo + turndowns.6mo| countyfips + year| 0 | 
                                     countyfips ,data = gr_corrected_scaled[gr_corrected_scaled$contra==1,])
  
  copart_fig <- as.data.frame(summary(gr_aug_corr_copart_stand)$coefficients[c("propertypercapl.6mo","disdecs.all.6mo","turndowns.6mo"),])
  copart_fig$var <- rownames(copart_fig)
  copart_fig$type <- "Co-Partisan"
  
  swing_fig <- as.data.frame(summary(gr_aug_corr_swing_stand)$coefficients[c("propertypercapl.6mo","disdecs.all.6mo","turndowns.6mo"),])
  swing_fig$var <- rownames(swing_fig)
  swing_fig$type <- "Swing" 
  
  contra_fig <- as.data.frame(summary(gr_aug_corr_contra_stand)$coefficients[c("propertypercapl.6mo","disdecs.all.6mo","turndowns.6mo"),])
  contra_fig$var <- rownames(contra_fig)
  contra_fig$type <- "Contra-Partisan"
  
  split_fig <- as.data.frame(rbind(copart_fig,swing_fig,contra_fig))
  split_fig$type <- factor(split_fig$type,levels=c("Co-Partisan","Swing","Contra-Partisan"))
  split_fig$var <- plyr::revalue(split_fig$var,c("propertypercapl.6mo"="Damage",
                                                 "disdecs.all.6mo"="Declarations",
                                                 "turndowns.6mo"="Turndowns"))
  split_fig$var <- factor(split_fig$var,levels=c("Damage","Declarations","Turndowns"))
  
  pdf("./results/gr_split_fig.pdf",width=8,height=5)
  print(
    ggplot(data=split_fig,aes(x=var,y=Estimate,ymin=Estimate-2*`Cluster s.e.`,ymax=Estimate+2*`Cluster s.e.`,group=type,colour=type))+
      geom_point(size=4,position=position_dodge(width=0.25))+
      geom_errorbar(size=2,width=0.25,position=position_dodge(width=0.25))+
      geom_hline(yintercept = 0,colour="indianred",size=2,linetype=2)+
      scale_colour_grey()+
      theme_minimal()+
      labs(colour = "Type of County")+
      xlab("Variable")+ylab("Standardized Coefficient")
  )
  dev.off()

  gr_split_sg <- stargazer(gr_aug_corr_copart,gr_aug_corr_swing,gr_aug_corr_contra,
                           dep.var.caption = c("Incumbent Party Two-Party Vote Share"),
                           dep.var.labels=c("Co-Partisan","Swing","Contra-Partisan"),
                           title="Disaster Damage and Incumbent Vote Share, 1972-2004: Split-Sample Estimates",
                           label="gr_split",
                           order=c(3,4,5,6,1,2),
                           covariate.labels = c("Damage",
                                                "Median Income",
                                                "Declarations, 6 Months",
                                                "Turndowns, 6 Months",
                                                "Incumbent Party Vote Share, Lag 1",
                                                "Incumbent Party Vote Share, Lag 2"
                           ),
                           keep.stat = c("n","rsq"),no.space = T,
                           star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                           notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                           table.layout ="-ld-t-as-n",
                           notes="\\parbox[t]{0.8\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       county-clustered standard errors in parentheses. All models include county and year fixed effects. 
                       $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(gr_split_sg, sep = '\n', file = "./results/gr_split.tex")

# models without FE or LDV
  
  gr_base_int_nofe <- felm(inc.twpct ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo*(copart+swing) +  medincK| year | 0 | countyfips ,data = gr_corrected)
  
  gr_aug_int_nofe <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+  disdecs.all.6mo + turndowns.6mo)*(copart+swing) +  medincK |  year| 0 | countyfips ,data = gr_corrected)
  
  gr_base_int_noldv <- felm(inc.twpct ~  propertypercapl.6mo*(copart+swing) +  medincK| countyfips + year | 0 | countyfips ,data = gr_corrected)
  
  gr_aug_int_noldv <- felm(inc.twpct ~  (propertypercapl.6mo+  disdecs.all.6mo + turndowns.6mo)*(copart+swing) +  medincK| countyfips + year| 0 | countyfips ,data = gr_corrected)
  
  gr_nofe_sg <- stargazer(gr_base_int_nofe,gr_aug_int_nofe,gr_base_int_noldv,gr_aug_int_noldv,
                          dep.var.labels = c("Incumbent Party Two-Party Vote Share"),
                          title="Impact of Disaster Damage on Incumbent Vote Share, 1972-2004: Alternative Models",
                          label="gr_nofe",
                          order=c(3,4,5,9,11,13,10,12,14,6,7,8,1,2),
                          covariate.labels = c("Damage",
                                               "Declarations, 6 Months",
                                               "Turndowns, 6 Months",
                                               "Co-Partisan X Damage",
                                               "Co-Partisan X Declarations",
                                               "Co-Partisan X Turndowns",
                                               "Swing X Damage",
                                               "Swing X Declarations",
                                               "Swing X Turndowns",
                                               "Co-Partisan",
                                               "Swing",
                                               "Median Income",
                                               "Incumbent Party Vote Share, Lag 1",
                                               "Incumbent Party Vote Share, Lag 2"
                          ),
                          add.lines = list(c("County Fixed Effects","","","\\checkmark","\\checkmark")),
                          keep.stat = c("n","rsq"),no.space = T,
                          star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                          notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                          table.layout ="-ld-#-t-as-n",font.size="footnotesize",
                          notes="\\parbox[t]{0.75\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       county-clustered standard errors in parentheses. All models include year fixed effects. 
                       $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(gr_nofe_sg, sep = '\n', file = "./results/gr_nofe.tex")
  

# split by party

  actual_inc_years <- c(1972,1976,1980,1984,1992,1996,2004)
  
  dem_inc_years <- c(1980,1996,2000)
  
  rep_inc_years <- c(1972,1976,1984,1988,1992,2004)

# run models in subsets of years with relevant parties

  gr_corrected$newout3 <- gr_corrected$inc.twpct
  gr_corrected$newout4 <- gr_corrected$inc.twpct
  
  gr_base_int_act <- felm(newout4 ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo*(copart+swing) +  medincK| countyfips + year | 0 | countyfips ,
                          data = gr_corrected[is.element(gr_corrected$year,actual_inc_years),])
  
  gr_aug_int_act <- felm(newout4 ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo)*(copart+swing) +  medincK| countyfips + year| 0 | countyfips ,
                         data = gr_corrected[is.element(gr_corrected$year,actual_inc_years),])
  
  gr_base_int_dem <- felm(inc.twpct ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo*(copart+swing) +  medincK| countyfips + year | 0 | countyfips ,
                          data = gr_corrected[is.element(gr_corrected$year,dem_inc_years),])
  
  gr_aug_int_dem <- felm(inc.twpct ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo)*(copart+swing) +  medincK| countyfips + year| 0 | countyfips ,
                         data = gr_corrected[is.element(gr_corrected$year,dem_inc_years),])
  
  gr_base_int_rep <- felm(newout3 ~ pres.vote.lag1  + pres.vote.lag2  + propertypercapl.6mo*(copart+swing) +  medincK| countyfips + year | 0 | countyfips ,
                          data = gr_corrected[is.element(gr_corrected$year,rep_inc_years),])
  
  gr_aug_int_rep <- felm(newout3 ~ pres.vote.lag1 + pres.vote.lag2 + (propertypercapl.6mo+disdecs.all.6mo + turndowns.6mo)*(copart+swing) +  medincK| countyfips + year| 0 | countyfips ,
                         data = gr_corrected[is.element(gr_corrected$year,rep_inc_years),])
  
  gr_party_sg <- stargazer(gr_base_int_act,gr_aug_int_act,gr_base_int_dem,gr_aug_int_dem,gr_base_int_rep,gr_aug_int_rep,
                           dep.var.caption = c("Incumbent Party Two-Party Vote Share"),
                           dep.var.labels=c("Incumbents","Democrats","Republicans"),
                           title="Impact of Disaster Damage on Incumbent Vote Share, 1972-2004: Split by Party",
                           label="gr_party",
                           order=c(3,6,7,9,11,13,10,12,14,4,5,8,1,2),
                           covariate.labels = c("Damage",
                                                "Declarations, 6 Months",
                                                "Turndowns, 6 Months",
                                                "Co-Partisan $\\times$ Damage",
                                                "Co-Partisan $\\times$ Declarations",
                                                "Co-Partisan $\\times$ Turndowns",
                                                "Swing $\\times$ Damage",
                                                "Swing $\\times$ Declarations",
                                                "Swing $\\times$ Turndowns",
                                                "Co-Partisan",
                                                "Swing",
                                                "Median Income",
                                                "Incumbent Party Vote Share, Lag 1",
                                                "Incumbent Party Vote Share, Lag 2"
                           ),
                           keep.stat = c("n","rsq"),no.space = T,font.size="footnotesize",
                           star.cutoffs = c(0.10,0.05),star.char = c("*","**"),table.placement="!ht",
                           notes.append = FALSE,notes.label = "",column.sep.width="0pt",
                           table.layout ="-ld-#-t-as-n",
                           notes="\\parbox[t]{0.9\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                       county-clustered standard errors in parentheses. All models include county and year fixed effects. 
                       $^{*}$p$<$0.10, $^{**}$p$<$0.05 (two-tailed test).}")
  
  cat(gr_party_sg, sep = '\n', file = "./results/gr_party.tex")
  

# write a stata file for the predicted values plot

  foreign::write.dta(gr_corrected,"gr_corrected_stata.dta")

# make data frames with mean values of everything (and just pick a year) but one is copartisan and another is not
# these values don't hugely matter because we're going to normalize the treatment effect with zero damage to 0

  newdata_contra <- as.data.frame(cbind(
    seq(0,max(gr_corrected$propertypercapl.6mo,na.rm=T),by=0.5),
    mean(gr_corrected$pres.vote.lag1,na.rm=T),
    mean(gr_corrected$pres.vote.lag2,na.rm=T),
    mean(gr_corrected$disdecs.all.6mo,na.rm=T),
    mean(gr_corrected$turndowns.6mo,na.rm=T),
    0,0,
    mean(gr_corrected$medincK,na.rm=T),
    1,0,0,0,0,0,0,0),
  )
  
  colnames(newdata_contra) <- c("propertypercapl.6mo",
                                "pres.vote.lag1",
                                "pres.vote.lag2",
                                "disdecs.all.6mo",
                                "turndowns.6mo",
                                "copart",
                                "swing",
                                "medincK",
                                "_Iyear_1976",
                                "_Iyear_1980",
                                "_Iyear_1984",
                                "_Iyear_1988",
                                "_Iyear_1992",
                                "_Iyear_1996",
                                "_Iyear_2000",
                                "_Iyear_2004")
  
  foreign::write.dta(newdata_contra,"newdata_contra.dta")
  
  newdata_copart <- newdata_contra
  newdata_copart$copart <- 1
  
  foreign::write.dta(newdata_copart,"newdata_copart.dta")

